// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

//! 
//!  Description: xyz mol filestream struct
//!  version 0.2
//!  typename T is 3d vector type and has to have  .x() .y() .z() members
//!
//!  

#ifndef XYZ_STREAM_H
#define XYZ_STREAM_H

#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <vector>
#include <set>
#include <map>
#include <iomanip>
#include <stdexcept>
#include <Eigen/Eigen>

namespace io
{

template<typename T>
struct xyz
{
public:
	struct record_type
	{
		short int mark;
		std::string label;
		T coord;
		T veloc;
		T force;
		bool operator < (const record_type & right) const
		{
			double tiny = 1e-8;
			if(coord.z() < right.coord.z()) return true;
			if(std::abs(coord.z() - right.coord.z()) <= tiny )
			{
				if(coord.y() < right.coord.y()) return true;
				if(std::abs(coord.y() - right.coord.y()) <= tiny )
					if(coord.x() < right.coord.x()) return true;
			}
			return false;
		}
	};
private:
	T a1_,a2_,a3_;
	std::vector<record_type> data_;
	std::map<std::string,double> double_var;
	std::map<std::string,int> int_var;
	std::map<std::string,std::string> string_var;
private:
	inline int index_converter(int index) const
	{
		int raw_index = -1;
		if(index<0) raw_index = data_.size()+index;
		if(index>0) raw_index = index-1;
		if((raw_index < 0) || (raw_index >= (int)data_.size()) )
			throw std::out_of_range("wrong xyz index");
		return raw_index;
	}
public:
	xyz() {}
	xyz(const T &a1,const T &a2,const T &a3,
		const std::vector<std::string> &label,
		const std::vector<T> &coord,const std::vector<T> &veloc,const std::vector<T> &force)
	{
		a1_ = a1;
		a2_ = a2;
		a3_ = a3;
		if( ( label.size() != coord.size()) || 
			( label.size() != veloc.size()) || 
			( label.size() != force.size())  )
			throw std::range_error("xyz set error coord veloc force");

		data_.resize(label.size());
		for(size_t i = 0; i < data_.size(); i++)
		{
			data_[i].label = label[i];
			data_[i].coord = coord[i];
			data_[i].veloc = veloc[i];
			data_[i].force = force[i];
		}
	}

	void set(const T &a1,const T &a2,const T &a3,
		const std::vector<std::string> &label,
		const std::vector<T> &coord,const std::vector<T> &veloc,const std::vector<T> &force)
	{
		a1_ = a1;
		a2_ = a2;
		a3_ = a3;
		if( ( label.size() != coord.size()) || 
			( label.size() != veloc.size()) || 
			( label.size() != force.size())  )
			throw std::range_error("xyz set error coord veloc force");

		data_.resize(label.size());
		for(size_t i = 0; i < data_.size(); i++)
		{
			data_[i].label = label[i];
			data_[i].coord = coord[i];
			data_[i].veloc = veloc[i];
			data_[i].force = force[i];
		}
	}

	void get_label(std::vector<std::string> & label)
	{	
		label.resize(data_.size());
		for(size_t i = 0; i < data_.size(); i++)
			label[i] = data_[i].label;
	}
	void get_coord(std::vector<T> & coord)
	{	
		coord.resize(data_.size());
		for(size_t i = 0; i < data_.size(); i++)
			coord[i] = data_[i].coord;
	}
	void get_veloc(std::vector<T> & veloc)
	{	
		veloc.resize(data_.size());
		for(size_t i = 0; i < data_.size(); i++)
			veloc[i] = data_[i].veloc;
	}
	void get_force(std::vector<T> & force)
	{	
		force.resize(data_.size());
		for(size_t i = 0; i < data_.size(); i++)
			force[i] = data_[i].force;
	}

	inline void swap(int i,int j)
	{
		i = index_converter(i);
		j = index_converter(j);
		std::swap(data_[i],data_[j]);
	}
	inline int length() const
	{
		return (int)data_.size();
	}
	inline record_type & operator[] (int index)
	{
		return data_[index_converter(index)];
	}
	inline const record_type & operator[] (int index) const
	{
		return data_[index_converter(index)];
	}
	
	double dist(int i,int j) const
	{
		i = index_converter(i);
		j = index_converter(j);
		double dx = data_[i].coord.x()-data_[j].coord.x();
		double dy = data_[i].coord.y()-data_[j].coord.y();
		double dz = data_[i].coord.z()-data_[j].coord.z();

		return sqrt(dx*dx+dy*dy+dz*dz);
	}

	inline void mark(const std::list<int> & index_set)
	{
		for(std::list<int>::const_iterator it = index_set.begin(); it != index_set.end(); it++)
			data_[index_converter(*it)].mark = 1;
	}
	
	inline std::vector<record_type> & data()
	{
		return data_;
	}
	inline const std::vector<record_type> & data() const
	{
		return data_;
	}

	inline T & a1()
	{
		return a1_;
	}
	inline const T & a1() const
	{
		return a1_;
	}
	inline T & a2()
	{
		return a2_;
	}
	inline const T & a2() const
	{
		return a2_;
	}
	inline T & a3()
	{
		return a3_;
	}
	inline const T & a3() const
	{
		return a3_;
	}

	inline double & var_double(const std::string & var_name)
	{
		return double_var[var_name];
	}
	inline std::string & var_string(const std::string & var_name)
	{
		return string_var[var_name];
	}
	inline int & var_int(const std::string & var_name)
	{
		return int_var[var_name];
	}

	inline const double & var_double(const std::string & var_name) const
	{
		return double_var[var_name];
	}
	inline const std::string & var_string(const std::string & var_name) const
	{
		return string_var[var_name];
	}
	inline const int & var_int(const std::string & var_name) const
	{
		return int_var[var_name];
	}

	inline void remove_var(const std::string & var_name)
	{
		double_var.erase(var_name);
		int_var.erase(var_name);
		string_var.erase(var_name);
	}

	void operator += (xyz & other)
	{
		for(typename std::vector<record_type>::iterator it = other.data_.begin(); it != other.data_.end(); it++)
			data_.push_back(*it);
	}
	
	void operator += (const record_type & add)
	{
		data_.push_back(add);
	}

	template<typename Ty>
	friend std::istream & operator >> (std::istream & is, xyz<Ty> & out);
	template<typename Ty>
	friend std::ostream & operator << (std::ostream & os, const xyz<Ty> & in);
	template<typename Ty>
	friend void print(std::ostream & os, const xyz<Ty> & in);
};




template<typename T>
std::istream & operator >> (std::istream & is, io::xyz<T> & out)
{
	out.data_.clear();
	std::string line;
	int length = 0;
	if(std::getline(is,line))
	{
		std::istringstream iss(line);
		iss >> length;
	}

	out.a1_.x() = 0.0;
	out.a1_.y() = 0.0;
	out.a1_.z() = 0.0;
	out.a2_.x() = 0.0;
	out.a2_.y() = 0.0;
	out.a2_.z() = 0.0;
	out.a3_.x() = 0.0;
	out.a3_.y() = 0.0;
	out.a3_.z() = 0.0;

	if(std::getline(is,line))
	{
		std::istringstream iss(line);
		double x,y,z;
		x = 0.0;
		y = 0.0;
		z = 0.0;
		iss >> x >> y >> z;
		out.a1_.x() = x;
		out.a1_.y() = y;
		out.a1_.z() = z;
		iss >> x >> y >> z;
		out.a2_.x() = x;
		out.a2_.y() = y;
		out.a2_.z() = z;
		iss >> x >> y >> z;
		out.a3_.x() = x;
		out.a3_.y() = y;
		out.a3_.z() = z;
		std::string str;
		while(iss >> str)
		{
			std::string::size_type pos = str.find_first_of("=", 0);
			str[pos]=' ';
			std::istringstream isst(str);
			std::string var_name;
			isst >> var_name;
			if(var_name[0] == '%')
			{
				int value;
				isst >> value;
				out.int_var[var_name] = value;
			} else if(var_name[0] == '$')
			{
				std::string value;
				isst >> value;
				out.string_var[var_name] = value;
			}
			else
			{
				double value;
				isst >> value;
				out.double_var[var_name] = value;
			}
		}
	}

	for(int i = 0; i < length; i++)
	if(std::getline(is,line))
	{
		std::istringstream iss(line);
		typename xyz<T>::record_type t;
		iss >> t.label;
		double x,y,z;
		iss >> x >> y >> z;
		t.coord.x() = x;
		t.coord.y() = y;
		t.coord.z() = z;
		double vx,vy,vz;
		vx = 0.0;
		vy = 0.0;
		vz = 0.0;
		iss >> vx >> vy >> vz;
		t.veloc.x() = vx;
		t.veloc.y() = vy;
		t.veloc.z() = vz;
		double fx,fy,fz;
		fx = 0.0;
		fy = 0.0;
		fz = 0.0;
		iss >> fx >> fy >> fz;
		t.force.x() = fx;
		t.force.y() = fy;
		t.force.z() = fz;
		t.mark = 0;
		out.data_.push_back(t);
	}

	if(length != out.data_.size()) throw std::length_error("Corrupt .xyz file.");
	return is;
};



template<typename T>
std::ostream & operator << (std::ostream & os, const io::xyz<T> & in)
{
	os << std::left << std::setw(5) << in.length();
	os << std::endl;

	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << in.a1_.x() << ' ';
	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << in.a1_.y() << ' ';
	os << std::left << std::setw(16) << std::setprecision(6) << std::fixed << in.a1_.z() << ' ';

	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << in.a2_.x() << ' ';
	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << in.a2_.y() << ' ';
	os << std::left << std::setw(16) << std::setprecision(6) << std::fixed << in.a2_.z() << ' ';

	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << in.a3_.x() << ' ';
	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << in.a3_.y() << ' ';
	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << in.a3_.z() << ' ';

	for(std::map<std::string,double>::const_iterator it = in.double_var.begin(); it != in.double_var.end(); it++)
		os << it->first << '=' << it->second << ' '; 
	for(std::map<std::string,int>::const_iterator it = in.int_var.begin(); it != in.int_var.end(); it++)
		os << it->first << '=' << it->second << ' '; 
	for(std::map<std::string,std::string>::const_iterator it = in.string_var.begin(); it != in.string_var.end(); it++)
		os << it->first << '=' << it->second << ' '; 
	os << std::endl;

	for(int i = 1; i <= in.length(); i++)
	{
		os << std::left << std::setw(5) << in[i].label;

		os << std::right << std::setw(12) << std::setprecision(6) << std::fixed << in[i].coord.x() << ' ';
		os << std::right << std::setw(12) << std::setprecision(6) << std::fixed << in[i].coord.y() << ' ';
		os << std::right << std::setw(12) << std::setprecision(6) << std::fixed << in[i].coord.z() << ' ';

		os << std::right << std::setw(16) << std::setprecision(6) << std::fixed << in[i].veloc.x() << ' ';
		os << std::right << std::setw(12) << std::setprecision(6) << std::fixed << in[i].veloc.y() << ' ';
		os << std::right << std::setw(12) << std::setprecision(6) << std::fixed << in[i].veloc.z() << ' ';

		os << std::right << std::setw(16) << std::setprecision(6) << std::fixed << in[i].force.x() << ' ';
		os << std::right << std::setw(12) << std::setprecision(6) << std::fixed << in[i].force.y() << ' ';
		os << std::right << std::setw(12) << std::setprecision(6) << std::fixed << in[i].force.z() << ' ';

		os << std::endl;
	}

	return os;
}




template<typename T>
void print(std::ostream & os, const io::xyz<T> & in)
{
	os << "unit cell vectors" << std::endl;
	os << std::left << std::setw(6) << std::setprecision(2) << std::fixed << in.a1_.x() << ' ';
	os << std::left << std::setw(6) << std::setprecision(2) << std::fixed << in.a1_.y() << ' ';
	os << std::left << std::setw(6) << std::setprecision(2) << std::fixed << in.a1_.z() << ' ';
	os << std::endl;

	os << std::left << std::setw(6) << std::setprecision(2) << std::fixed << in.a2_.x() << ' ';
	os << std::left << std::setw(6) << std::setprecision(2) << std::fixed << in.a2_.y() << ' ';
	os << std::left << std::setw(6) << std::setprecision(2) << std::fixed << in.a2_.z() << ' ';
	os << std::endl;

	os << std::left << std::setw(6) << std::setprecision(2) << std::fixed << in.a3_.x() << ' ';
	os << std::left << std::setw(6) << std::setprecision(2) << std::fixed << in.a3_.y() << ' ';
	os << std::left << std::setw(6) << std::setprecision(2) << std::fixed << in.a3_.z() << ' ';
	os << std::endl;

	for(std::map<std::string,double>::const_iterator it = in.double_var.begin(); it != in.double_var.end(); it++)
		os << it->first << '=' << it->second  << std::endl;; 
	for(std::map<std::string,int>::const_iterator it = in.int_var.begin(); it != in.int_var.end(); it++)
		os << it->first << '=' << it->second << std::endl;; 
	for(std::map<std::string,std::string>::const_iterator it = in.string_var.begin(); it != in.string_var.end(); it++)
		os << it->first << '=' << it->second << std::endl;; 

	for(int i = 1; i <= (int)in.length(); i++)
	if(in[i].mark == 1)
	{
		os << std::left << std::setw(4) << i;
		os << std::left << std::setw(5) << i - 1 - (int)in.length();
		os << std::left << std::setw(4) << in[i].label;

		os << std::right << std::setw(6) << std::setprecision(2) << std::fixed << in[i].coord.x() << ' ';
		os << std::right << std::setw(6) << std::setprecision(2) << std::fixed << in[i].coord.y() << ' ';
		os << std::right << std::setw(6) << std::setprecision(2) << std::fixed << in[i].coord.z() << ' ';

		os << std::right << std::setw(8) << std::setprecision(2) << std::fixed << in[i].veloc.x() << ' ';
		os << std::right << std::setw(6) << std::setprecision(2) << std::fixed << in[i].veloc.y() << ' ';
		os << std::right << std::setw(6) << std::setprecision(2) << std::fixed << in[i].veloc.z() << ' ';

		os << std::right << std::setw(8) << std::setprecision(2) << std::fixed << in[i].force.x() << ' ';
		os << std::right << std::setw(6) << std::setprecision(2) << std::fixed << in[i].force.y() << ' ';
		os << std::right << std::setw(6) << std::setprecision(2) << std::fixed << in[i].force.z() << ' ';

		os << std::endl;
	}
}


}
/*
//! coord expression class 
//! coord expression is itmes separated by delimiters like a,b,c
//! note that coord expression cannot contain whitespace
class coord_expression
{
public:
private:
	Eigen::Vector3d v;
	const std::map<int,Eigen::Vector3d> & index_map;
public:
	//! constructor
	//! @param d list delimiter (,)
	coord_expression(const std::map<int,Eigen::Vector3d> & index_map) : index_map(index_map) {}
	//! * operator returns with the list reference
	//! @return list_type const reference
	const Eigen::Vector3d & operator* () const { return v; }
	//! [] operator can give const access to the list items by their index [0..length-1]
	const T & operator[] (size_t index) const { return l[index]; }
	
	//! istream operator
	friend std::istream & operator >> (std::istream & is,coord_expression & expr);
};

//! istream operator for coord_expression
std::istream & operator >> (std::istream & is,coord_expression & expr)
{
	std::string m;
	if ( !(is >> m) ) return is;
	expr.l.clear();
	std::replace(m.begin(),m.end(),expr.d,' ');
	std::istringstream isd(m);
	T item;
	while(isd >> item) expr.l.push_back(item);
	return is;
}
*/

#endif // XYZ_STREAM_H
